The main data used in this tutorial and in the lecture are about the geolocalisation of french restaurants in Paris. Data is extracted from an official register called SIRENE (Computer system for the business and establishment register) managed by the French National Institute of Statistics and Economic Studies (Insee) and geolocated by Etalab (French task force for Open Data). This register records the civil status of all companies and their establishments (including restaurants). SIRENE has the advantages of being rigorous and exhaustive on the French territory.
sf::st_read().
Reading layer `iris_75' from data source `/home/tim/Documents/prz/satRday/exercises/data/iris_75.shp' using driver `ESRI Shapefile'
Simple feature collection with 992 features and 4 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 643075.6 ymin: 6857477 xmax: 661086.2 ymax: 6867081
epsg (SRID): NA
proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
plot(iris_75). What do you notice ?
sf::st_geometry() function? What solution do you propose then?
sf::st_geometry() makes it possible to isolate the information contained in the ‘geometry’ column of the sf object. Using it, we put aside other variables (here CODE_IRIS, P14_POP, AREA and CODE_COM).
st_read() and sf::st_geometry().
Reading layer `sir_75' from data source `/home/tim/Documents/prz/satRday/exercises/data/sir_75.shp' using driver `ESRI Shapefile'
Simple feature collection with 23348 features and 7 fields
geometry type: POINT
dimension: XY
bbox: xmin: 643502.1 ymin: 6857645 xmax: 659766.4 ymax: 6867041
epsg (SRID): NA
proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
plot(st_geometry(iris_75), bg = "cornsilk", col = "lightblue",
border = "white", lwd = .5)
plot(st_geometry(sir_75), col = "red", pch = 20, cex = .2, add=TRUE)
title("Restaurants in Paris")sf::st_intersects() and sapply().
wwww
inter <- st_intersects(x = iris_75, y = sir_75)
iris_75$RESTAU <- sapply(inter, length)
head(iris_75)Simple feature collection with 6 features and 5 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 650181.1 ymin: 6861761 xmax: 652179 ymax: 6863138
epsg (SRID): NA
proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
CODE_IRIS P14_POP AREA CODE_COM geometry RESTAU
1 751010101 974 64006.07 75101 MULTIPOLYGON (((652096.5 68... 44
2 751010102 167 62671.06 75101 MULTIPOLYGON (((652052.7 68... 6
3 751010103 246 47653.41 75101 MULTIPOLYGON (((651869 6862... 16
4 751010104 3 222656.99 75101 MULTIPOLYGON (((651639.9 68... 10
5 751010105 0 243255.23 75101 MULTIPOLYGON (((650369.8 68... 0
6 751010199 0 234277.01 75101 MULTIPOLYGON (((652096.5 68... 0
dplyr package: select, group_by et summarize. These functions also work with sf objects.
We would like here to design EPCI maps that combine the number of restaurants and the number of restaurants per 10,000 inhabitants. The EPCI correspond to the level of french intercommunalites.
Data preparation:
getBreaks et carto.pal functions of the cartography package. For the creation of the typo variable, you can use the cut function and apply the parameters digit.lab = 2 and include.lowest = TRUE.
library(sf)
library(cartography)
library(dplyr)
# Import data
fra <- st_read("data/fra.shp", quiet = TRUE)
epci_31 <- readRDS("data/epci_31.rds")
# Create the variable
epci_31$nb_rest_10000inhab <- 10000 * epci_31$nb_of_rest / epci_31$P14_POP
# Define breaks
bks <- getBreaks(v = epci_31$nb_rest_10000inhab, method = "quantile", nclass = 4)
# Define color palette
cols <- carto.pal("orange.pal", length(bks)-1)
# Create a "typo"" variable
epci_31 <- epci_31 %>%
mutate(typo = cut(nb_rest_10000inhab,breaks = bks, dig.lab = 2,
include.lowest = TRUE))cartography package, make the following map which contains in a choropleth layer the variable nb_rest_10000inhab and in a proportional circle layer the variable nb_of_rest. The you can try to do the same map using the ggplot2 and tmaps packages.
With cartography:
# Define plot margins
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
# Find EPCI bounding box
bb <- st_bbox(epci_31)
# Plot France using EPCI boundingbox
plot(st_geometry(fra), col="ivory", border = "ivory3",
xlim = bb[c(1, 3)], ylim = bb[c(2, 4)])
# Plot the choropleth layer
choroLayer(epci_31, var = "nb_rest_10000inhab",
breaks = bks, col = cols, border = "grey80", lwd = 0.5,
legend.pos = "topleft",add = TRUE,
legend.title.txt = "Number of restaurants\nfor 10,000 inhabitants")
# Plot proportionnal symbols
propSymbolsLayer(epci_31, var="nb_of_rest", col="#440170",border=NA,
legend.pos="left", inches=0.4, add = TRUE,
legend.title.txt = "Number of restaurants")
# Add a layout layer
layoutLayer(title = "Restaurants", sources = "Insee, 2018",
author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred",
coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
# Add a north (south) arrow
north(pos = "topright", south = TRUE)With ggplot2:
library(ggplot2)
map_ggplot <- ggplot() +
geom_sf(data = fra, colour = "ivory3",
fill = "ivory") +
geom_sf(data = epci_31, aes(fill = typo), colour = "grey80") +
scale_fill_manual(name = "Number of restaurants\nfor 10,000 inhabitants",
values = cols, na.value = "#303030")+
geom_sf(data = epci_31 %>% st_centroid(),
aes(size= nb_of_rest), color = "#440154CC", show.legend = 'point')+
scale_size(name = "Number of restaurants",
breaks = c(1, 500, 3200),
range = c(0,18))+
coord_sf(crs = 2154, datum = NA,
xlim = st_bbox(epci_31)[c(1,3)],
ylim = st_bbox(epci_31)[c(2,4)]
) +
theme_minimal() +
theme(panel.background = element_rect(fill = "azure",color=NA)) +
labs(
title = "Restaurants",
caption = "Insee, 2018\nKim & Tim, 2018"
)
plot(map_ggplot)With tmap:
TODO Timothée
TODO Timothée
mapview package. Try using different parameters to customize your map.
map.types, col.regions, label, color, legend, layer.name, homebutton, lwd … parameters of the mapview function.
library(mapview)
library(sf)
sir_31 <- readRDS("data/sir_31.rds")
mapview(sir_31, map.types = "OpenStreetMap",
col.regions = "#940000",
label = paste(sir_31$L2_NORMALISEE, sir_31$NOMEN_LONG, sep = " - "),
color = "white", legend = TRUE, layer.name = "Restaurants in SIRENE",
homebutton = FALSE, lwd = 0.5)